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Abstract 

The Taylor Interpolation through FFT (TFFFT) algorithm for the computation of the elec- 
tromagnetic wave propagation in the quasi-planar geometry within the half-space is proposed in 
this article. There are two types of TFFFT algorithm, i.e., the spatial TFFFT and the spectral 
TFFFT. The former works in the spatial domain and the latter works in the spectral domain. It 
has been shown that the optimized computational complexity is the same for both types of TFFFT 
algorithm, which is AT° pt Af° pt O(Nlog 2 N) for an N = Af x x J\f y computational grid, where 7V° pt is 
the optimized number of slicing reference planes and N° pt is the optimized order of Taylor series. 
Detailed analysis shows that Af° pt is closely related to the algorithm's computational accuracy 7 TI , 
which is given as Af£ pt ~ — hi7 TI and the optimized spatial slicing spacing between two adjacent 
spatial reference planes 5° pt only depends on the characteristic wavelength A c of the electromagnetic 
wave, which is given as 5° pt ~ jj^c- The planar TFFFT algorithm allows a large sampling spacing 
required by the sampling theorem. What's more, the algorithm is free of singularities and it works 
particularly well for the narrow-band beam and the quasi-planar geometry. 
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I. INTRODUCTION 



The computation of electromagnetic wave propagation using the direct integration 
method is not efficient for the large-scale computation because the direct integration method 
has a daunting computational complexity of O {N 2 ) for an iV = Af x x Af y computational 



grid, e.g., in the beam-shaping mirror system design 



or the Quasi-Optical (QO) gyrotron 
Fortunately, when the computa- 



application, days of computation is required [l|, 3 Q |j 0] 
tional geometry is a plane, the FFT has been shown to be efficient in the electromagnetic 
wave computation 0,0,0], which has a computational complexity of 0(iVlog 2 N) and a low 
sampling rate only limited by the Nyquist rate. For the quasi-planar geometry, it will be 
shown in this article that the FFT can still be used with the help of the Taylor Interpolation 
(TI) technique. 

The rest of this article is organized as follows. Section |H] gives the 2-Dimensional (2D) 
Fourier spectrum of the electromagnetic wave in its closed-form expression. Section IIHI 
presents the optimized spatial and spectral types of TI-FFT algorithm. In Section IIV[ 
one numerical example is used to show the performance of the planar TI-FFT algorithm. 
Section |V| discusses the advantages and problems of the planar TI-FFT algorithm; some 
helpful suggestions are given. Finally, Section IVT1 summarizes the planar TI-FFT algorithm. 
The scheme used to illustrate the planar TI-FFT algorithm is shown in Fig. ^ and the time 
dependence e J ' w * has been assumed in this article. 

II. ELECTROMAGNETIC WAVE IN THE SPECTRAL DOMAIN 

In this section, the 2D Fourier spectrum and far-field of the electromagnetic wave for the 
radiation integral are shown to be closely related to each other. 

A. The radiation integral 

For given electric and magnetic surface currents (J s ,J ms ), the radiating electric field E 
can be obtained under the Lorenz condition 0, , which is given as 



(jJ6 JJS 



( 



k 2 J s (r>)G(R) 



J s (r') ■ V V'G(R) -jue3 ms (r') x V'G(R) 



V 



dV, 
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FIG. 1: Electromagnetic wave propagation and scattering: the computation of the electromagnetic 
wave propagation (the incident file E l ) onto the PEC surface S is implemented through the spatial 
TI-FFT and the computation of the scattered electromagnetic field E s from the PEC surface S 
is implemented through the spectral TI-FFTs and the inverse Fourier transform. 5z is the spatial 
slicing spacing in the spatial TI-FFT. 

where, V' is the gradient operator on the source coordinate r' and the scalar Green's function 
is given as 



G(R) 



47r|R| ' 



R = r - r. 



(2) 



B. The 2D Fourier spectrum of the scalar Green's function 



Now apply the 2D Fourier transform on the scalar Green's function G(R) in (J2J), 





where k z and the 2D Fourier transform has been defined as 




ki + e y <e 
kl+ki>e 



(4) 
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y=-oo 



e jkyV dy^ dx, 
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From ©, 
Q[k X i ky, r ) 



I /*oo 
e jk x x' e jk y y' j ^ e jk x (x-x') x 
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e — jfc|R| 

47r|R| 
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dx, (6) 



FT, 



G(R) 



2tt 



jk x x' ^jk y y' 
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and w = z — z', 


(jHJ) becomes, 
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J u= — oo 


4tt|R| 



> du (7) 



In the cylindrical coordinate, 



|R| = v /( r± )2 +w 2 



(8) 



where r± — u + v and the following relation can be obtained from 



dr 



JRj 



d|R| 



(9) 



Now, express in the cylindrical coordinate with the help of 
Q(k x , k y , r ) 



1 



47T 



1 /" 27r 



where t/> = arctan 



R|=|to| 



and d> = arctan 



27T Jd>=0 



-jk±r± cos(i/>- 



d|R| (10) 



. The integration over is the Bessel function of 



the first kind of order and (jlUj) reduces to 



Q{k x , ky, r ) 



gjkxX ^jk y y 



n\=\w\ 



e -jk\R\ x j 



jfei a/IRP-w; 2 



d|R| 



_ e 3"feaa!' e jfc w i/ / e -ifc Jt |«-« / | 



47rA; ? 



(11) 



Because only half-space 2 > 2' is of interest, only the 2D Fourier spectrum for half-space 



z > z' will be considered in the rest of this article, which is obtained from (fTTl) as 



Q > {k x , ky, r ) 



Airk 



(12) 
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C. 2D Fourier spectra of Green's function related expressions 



The 2D Fourier spectra of the derivatives (order n 
obtained from the property of the Fourier transform 



of the scalar Green's function can be 



a 



d^G(R) 



/ -jk T ) n g > (k x ,k y ,r f ), r = x,y,z. 



(13) 



Particularly, for the first-order and second-order derivatives, 



(14) 



T^^ e * T=x > y >*- (15) 

Similarly, the 2D Fourier spectra of the following expressions can be obtained for half- 
space z > z', 



FT, 



VG(R) 



-jk.g > (k x ,k y ,r'). 



(16) 



FT 



V 2 G(R) 



-k Q > {k x , ky, r ). 



(17) 



FT, 



VVG(R) 



kkQ > [k x , ky, r ) 



(18) 



FT, 



GJR) 



Q > [kx-, ky, r ) 
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(19) 



FT, 



Gm.{ R) 



kx I 



(20) 



where the dyadic Green's functions of the electric type (G e ) and the magnetic type (G Tl 
are given as 



G e (R) = (j+^Vv) G(R). 



(21) 



G m (R) = VG(R) x I. 



(22) 
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D. The far-fields 



In the far-field limit, R ~ r — > oo, 

G(r) = ( — 1 ^ r , R = r in the far-field limit. (23) 
47r|r| 

Similarly, the first-order derivative of the Green's function in the far-field limit can be 
obtained as 



dG(r) t ( 1 \ e' jk ^ ( r , \ e^' fc l r l , e~ jk \ r 

— ^r~^ = i — r ~~ Jk 2 — r~r — ~J tt^ — r~r = ~jK — tt 

ot r V Irl / 47r r V r / Air r 47T r 



where only jij term is kept and the other terms (t^, nj, • • •) are ignored. In derivation of 
(|24jl . the following relation has been used in the far-field limit, 

— = pr, r = x,y,z. (25) 

K |r| 

Following the similar procedure given in (|24|). the far- fields of derivatives (order n) of the 
scalar Green's function are obtained as 

d^Glr) e -J k \ r \ 

~d^r = HK)n ^\' r = x,y,z. (26) 

It is not difficult to see that the far-fields and the 2D Fourier spectra are closely related 
to each other. 



E. The 2D Fourier spectra of the 3D spatial convolutions 

It is not difficult to see that, the radiation integral in (Q) can be expressed as the sum of the 
3D spatial convolutions of some source terms with the Green's function related expressions. 
For simplicity, let's consider the 3D spatial convolution of an arbitrary source term s with 
the scalar Green's function G, 



3D 

sir) 



(g) G(r) = s(r')G(R) dS', (27) 

Now, apply the 2D Fourier transform on (|2~7jl and express the scalar Green's function 
G(R) in the spectral domain, 



S(k x ,k y ) = FT 2 



s(r) (g) G(r) 



(28) 



2,71 Jx=—oo Jy=—oa 



s r x 



2,71 Jk' x =—oo Jk' y =—oo 



-jk' x (x-x') -jk' y {y-y') 



x - — dk~dk. 



dS' > dxdy, 



First, do the integral over (x,y), (|58|l reduces to 



S(k x ,k y ) 



sir 

S I ' Jk' 



oo roc 



-oo Jk' y =— oo 



gjkxx'gjkyy' 



(29) 



j'g jk'z( z z ') 

x — — S(k' x — k x )8{k' y — k y )dk' x dk' y 



AirkL 



dS' 



Next, do the integral over (k x , k' ) and (J2TIJ) reduces to 



FT, 



s(r)0G(r) =L s(r) ^(fc^^O), 







where L in ()3())) is the radiation vector |l£J| for source term s, which is defined as 



(30) 







s(r) 


"//. 







J J s{r')e* r 'dS', 



(31) 



It is not difficult to see that the radiation vector L in (13 1|) reduces to the regular 2D 
Fourier spectrum when surface S is a plane located at z' = 0. 



L s(r) 



2ttFT, 



z'=0 



s(x,y) 



(32) 



where the dummy primed (x',y f ) have been replaced with (x,y). Substitute (j31j) into ([30)1 . 
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S(k x ,kj 



jk z z 



4vrk 7 



s r 



£(k x ,k y ,0)L s(r) 



(33) 



V 



It is not difficult to see that the radiation vector L in (}3~T]) is closely related to the far-field 
by letting r — > oo, 



s(r) ® G(r) 



-jfc|r| 



/ 



47T |r I 



s r 



(34) 



V 



From (|34p. if ^j^r- can be ignored, the radiation vector L can be considered as the far- field 
pattern, which means that when the far-field is obtained, the radiation vector and the 2D 
Fourier spectrum of the 3D convolution are also obtained, from (J34)) and (J33)) respectively. 



F. 2D Fourier spectrum of the radiation integral 

From (|33|) and (JTJ), the 2D Fourier spectrum of the radiation integral (denoted as J 7 ) is 
obtained, which is 



T= -^P(k x ,k y ,0) 

UJ6 



k 2 L J(r) -k £ 



r=x,y,z 



k r L J T (r) 



JmW X k 



(35) 



G. Electromagnetic field on a plane 



After the 2D Fourier spectrum T has been obtained, the electric field E can be expressed 



in the PWS form 



in. 



12j . which is given as, 



Mr) = TFT, 



J- (k x ,k y )e-^ z 



(36) 



where Tq and the 2D Inverse Fourier Transform have been defined as 



.Fo(k x , k y ) = J" 0x x + T Qy y + T^i = J"(k x , k y )ej 



k z z 



s 



Oz 



ky 
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III. THE PLANAR TI-FFT ALGORITHM 



In this section, the optimized spatial and spectral TI-FFTs are presented. It will be 
shown that both of them have the same computational complexity for the same quasi-planar 
surface. 



A. The spatial TI-FFT algorithm 

It has been shown in (jHUj) that the electric field E on a plane can be evaluated through 
the 2D inverse Fourier transform. For a quasi-planar surface, the TI technique can be used, 
which leads to the spatial TI-FFT algorithm (where the quasi-planar surface is sliced into 
many small spatial subdomains, as shown Fig. [TJ. 

Rewrite the electric field E in (J55J) as follows, 



Efr) = e-^IFT, 



.Fo(kx,ky) 



JAk z Az 



■F (k x , k y ) = ^(kx, k y )eJ 



jZ\k z Zrrrir] 



(37) 



where z mln denotes the minimum value of z. Now, express e i Ak * Az into a Taylor series on 
the spatial reference plane located at z = z r , 



JAk z Az _ J&kz&zr 



Mo r 

E 

n=0 



— (jAk z ) n (z- z, 



m 



(38) 



where Az r = z r — z min and Af Q is the order of Taylor series. 

Substitute (|3*Hj) into (|37|h the spatial TI-FFT algorithm for the electric field E is obtained, 



EM 



Mo 

-jkz 

n=0 



1 

n! 



J [z — z r 



IFT, 



(39) 
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The number of spatial reference planes M T required in the computation depends on the 



spatial slicing spacing I 5 Z = max 



Zr+l 



and the characteristic surface vari- 



ation Az c (within which the electromagnetic field is of interest): J\f r oc Az c /5 Z . The readers 
should note that the actual maximum interpolation distance is y, which is located at the 
middle of two adjacent spatial reference planes, but 5 Z is used in this article to simplify the 
notation. Apparently, to achieve the desired computational accuracy (denoted as 7 TI ), the 
choice of the spatial slicing spacing 5 Z between two adjacent spatial reference planes depends 
on Ak zc , which is defined as 



k 



k 2 -k 2 



ka, 



1 



\ 



(40) 



where, fcj_ iC is the characteristic bandwidth (beyond which the 2D Fourier spectrum T is 
negligible) and is defined on x-y plane. It is clear that the smaller the bandwidth the 
larger the S z could be, which also means a smaller J\f r . So, a narrow-band beam and a small 
surface variation Az c (quasi-planar geometry) are in favor of the planar TI-FFT algorithm. 

In view of the importance of the spatial slicing spacing 5 Z , it is helpful to define the 
characteristic wave length A c for a narrow-band beam. From (J4*U|) . 



2tt 



2tt 



A 



k-Jk 2 -k 2 ±c «' 



(41) 



For a narrow-band beam (k± >c -C fc), 



A c ~2 



A. 



(42) 



Fig. |21 plots the exact value in (|41|) and approximation in (}4*2*|) of the characteristic 
wavelength A c for different characteristic bandwidth k± yC , from which it can be seen that 
the maximum deviation of the approximation from the exact value is 1A, which occurs at 
k± jC = k. 

It can be seen from and that, for the given computational accuracy j TI , the 
spatial slicing spacing 5 Z should satisfy the following relation, 



7ti ~ O 



Ak z c 5 Z 



Mo 



(43) 
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k, /k 

FIG. 2: The plots of the characteristic wavelength A c for different k ZtC . The exact value (line) 
is given in (|41j) and the approximation (dots) is given in Q42j) , The plots show that A c S> A for 
a narrow-band beam. The maximum deviation of the approximation from the exact value is 1A, 
which occurs at k Z)C = (&_i_ jC = k). 



5 ~ — ( — ) ^ = — ( — ) ™ (44) 
ka I 7 TI J 2ira 1 7ti / 



For a narrow-band beam (k± >c <C k). 



1 / k V ( 1 \ & 



^UdUr 1 A (45) 



/TI 



Now consider a quasi-planar surface with a characteristic surface variation of Az c = J\f z \, 
from (|44jl the number of spatial reference planes A/" r is given as 

i 

K = ~ 2™ — A4, (46) 
For a narrow-band beam (fc_i_ )C <C fc), 
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M r ~ 7T 



'jfe 



1 \S 



A4 



(47) 



k J \ 7 TI y 

The number of FFT operations N FFT and the computational complexity CPU are obtained 



as 



/ 1 \a7^ 

N FFT =M xM r = 2vra — M Q M Z 

\7ti J 



(48) 



CPU = N FFT O 
For a narrow-band beam (fcj_ )C -C k) 





( 1 \^ 




Nlog 2 N 


= 2vra f — j MoMz O 


N\og 2 N 









(49) 



k 



71 



1 \ ^° 



7t 



MoM z - 



(50) 



CPU 



A- 



~ 7T 



-L,c 



1 \K 



iVlog 2 iV 



(51) 



V k J \ 7 TI J 

For a narrow-band beam, the computational complexity CPU has a square law depen- 
dence on the characteristic bandwidth fcj_ )C of the electromagnetic wave and have a linear 
dependence on the surface variation (Az c = M Z X). The computational complexity CPU also 
has an inverse M th -root dependence on the computational accuracy 7 TI . So the character- 
istic bandwidth fcj_ )C has the most significant effect on the computational complexity of the 
planar TI-FFT algorithm. 

It can be seen from (J48j) or (|49|) that the optimized number of Taylor series M° pt can be 
obtained through finding the minimum value of N FFT in (|48|) or CPU in (j49J) by assuming 
that M is a continuous variable, 



dN F 



dM Q 







In [A/"o] -In [ 7l 



MT ~ round 



In 
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7ti 



dM 



round 
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FIG. 3: The linear dependence of the optimized order of Taylor series A/"o pt on the computational 
accuracy 7 TI (dB). It can be seen that N° pt = (2, 5, 7, 9) for 7 TI = (-20, -40, -60, -80) dB respec- 
tively. 

where "round" means to round the value to its nearest integer (actually, to achieve a higher 
computational accuracy j T1 , the upper-bound could be used but the computational com- 
plexity CPU is a little higher). Fig. El also shows the linear dependence of the optimized 
order of Taylor series Af° pt on the computational accuracy 7 TI (dB), which has been shown 
in (|53j). 

The optimized spatial slicing spacing 5° pt can be obtained from (}4"4"|) and (p)3*|) . which is 



where e ~ 2.718 is the natural logarithmic base. It is interesting to note that the optimized 
spatial slicing spacing 5 Z doesn't depend on the computational accuracy 7 TI and strongly 




(54) 



For a narrow-band beam (k±^ c k), 




(55) 
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FIG. 4: Plots of the exact value (line) of <5°. pt given in ([54)1 and approximation (dots) given in (155)) 
for different characteristic bandwidth k± tC . 

depends on the characteristic bandwidth k± )C (inverse square law). Fig. 0] shows <5° pt for 
different characteristic bandwidth k± tC , from which it can be seen that 5° z pt > 0.5A for 
k ZjC > 0.9k (jfe ±>c < 0.436&). 

The optimized number of spatial reference planes 7V° pt is given as 



A/l op 



AZr. 



2neaM z ~ 17aJ\f z , 



(56) 



For a narrow-band beam (/c_l jC <C k), 



A/I op 



'k 



~ lie 



±,c 
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(57) 



Substitute (J53|) into (J48|) and (J49|) . the optimized number of FFT operations iV° pt T and 
the optimized computational complexity CPU° pt can also be obtained, 



iVppIp = 2vrea In 



1 

7ti 



(58) 
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For a narrow-band beam (k± >c -C k), 



( k rh 


' i " 




7ti 



(59) 



The optimized computational complexity CPU opt is given as 



cpu opt ~ n; p ; t o 



N log 2 N 



(60) 



B. The spectral TI-FFT algorithm 

It has been shown in that the computation of the 2D Fourier spectrum T is equivalent 
to evaluate the radiation vector L. For the quasi-planar geometry, the FFT can still be used 
with the help of the TI technique, which leads to the spectral TI-FFT algorithm (where the 
spherical spectral surface is sliced into many small spectral subdomains, as shown Fig. EJ. 

From (|31|). the radiation vector L can be rewritten as 



f(r) 

V / 



27re Az ™FT 9 



f(r)e jk * Az 



(61) 



where f(r) = ^ and n is the normal to surface S. Now the Taylor expansion of L in ()61|) 
over k z is given as 

/ 



Mo 

f (r) | = 2ne jk * z ^ ^ < 

n=0 



, ( j [^z k z 



FT, 



f(r) Az 



(62) 



where k z ^ r denotes the spectral reference plane. For the given computational accuracy 7 TI , the 



spectral slicing spacing ( Sk z = max 
relation, 



k z k 7 _ 



k z , r +i — k z ^ should satisfy the following 



7ti 



(63) 
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(64) 
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k x + ky + k z — k 




FIG. 5: The spectral domain division for the spectral and TI-FFT: Only k z > half sphere surface 
is used for half-space z > z' computation in this article. k z ,r and kz,r+i denote the r th and (r + l) th 
spectral reference planes respectively. 5k z is the spectral slicing spacing. 

The number of spectral reference planes M r is given as 

K = =^~2ira — M z , (65) 
The number of FFT operations is given as 

/ i \ik 

A/" FPT ~ 2rta f — j A/"oA/" 2 . (66) 

It is obvious that 7V r in (JHSI) and J\f FFT in (jr]H|) are the same as those given in PrJj) and 
(I48|) . which also means that the spatial and spectral TI-FFTs have the same optimized 
computational complexity. 
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FIG. 6: The x-component (magnitude) of the scattered output field. 
IV. COMPUTATIONAL RESULTS 

To show the efficiency of the planar TI-FFT algorithm, the direct integration of the radi- 
ation integral in (JIJ has been used to make comparison with the planar TI-FFT algorithm. 
The numerical example used for such purpose is a 110 GHz Fundamental Gaussian Beam 
(FGB) scattered by a PEC quasi-planar surface with a sin wave perturbation. The 110 GHz 
FGB has a wavelength of A ~ 2.7 mm. 



A. The numerical results 

The incident 110 GHz FGB propagates at — z direction and has a beam waist radius of 
w — 1 cm. The quasi-planar PEC surface with a sine wave perturbation is described as 



z (x, y) = -2.5A + 0.5A cos (^tt— J cos ( 27r T^J • ( 67 ) 
In the numerical implementation of the planar TI-FFT algorithm, the computational 
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FIG. 7: The y-component (magnitude) of the scattered output field. 

accuracy j T1 = 0.0001 (—80 dB) has been used and the following optimized quantities are 
obtained from (Jo"3|) - (joTH) . 



N° pt ~ 9, 5T ~ 0.6A, AC° 



0.6 



NZl ~ 18, CPU 



opt 



180 



N\og 2 N 



.(68) 



where the quasi-planar surface described in (JBTj) has a characteristic surface variation Az c ~ 
1A. 

The scattered output field E s are evaluated on plane z = (where the incident 110 GHz 
FGB starts to propagate). Fig. |Hl Fig. [Hand Fig. Elshow the magnitude patterns of x-, y-, 
and z-components of the scattered output field E s . The comparison of the result obtained 
from the planar TI-FFT algorithm and that from the direct integration method is given in 
Fig. for both the magnitudes and the real parts, which shows that the planar TI-FFT 
algorithm has the desired —80 dB computational accuracy. 
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FIG. 8: The z-component (magnitude) of the scattered output field. 



B. The CPU time and the accuracy 



The CPU time t TI for the planar TI-FFT algorithm and t DI for the direct integration 
method have been summarized in Table El together with the coupling coefficient defined as 



C T = 



JJ E* [E* ]* dxdy 



(69) 



2 = 0, 



lH\K h r\ 2 dxdy^H\K,,r\ 2 dxdy 

where, E* r and E* : (r = x, y, z) denote the scattered output field components obtained 
from the planar TI-FFT algorithm and the direct integration method respectively. From 
TABLE HI it can be seen that, even though at a large sampling spacing 5 = 0.46A (M x = 
M y ~ 128), the coupling coefficients are still well above 90.00%. At this sampling rate, the 
direct integration method using Simpson's 1/3 rule is not accurate enough pi Q- [J- Also 
note that the coupling coefficients C T (r = x, y, z) reach their maximum values of 99.99% at 
A4 = My ~ 256 (8 X = 5 y = 0.23A), after which the accuracies remain constant and thus the 
Nyquist rate can be estimated roughly as M Nyquist ~ 256. The reason for this phenomenon 
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FIG. 9: The comparison of the scattered output field on plane z = 0, across the maximum value 
point of |E*| and at x direction: a) is for E*; and b) is for E s z ; solid lines (TI-FFT) and circles 
(direct integration method) are magnitudes; dashed lines (TFFFT) and dots (direct integration 
method) are real parts. 

is, that after the sampling rate increases above the Nyquist rate, further increasing the 
sampling rate will not give more information or computational accuracy. 

The CPU time for the planar TI-FFT algorithm t TI and for the direct integration method 



TABLE I: CPU time (t TI , t DI ) and coupling coefficient C T 
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FIG. 10: The CPU time (t TI , t DI ) comparison: bars in gray color are for the planar TI-FFT 
algorithm and bars in black color are for the direct integration method. Note that the CPU time 
is in logarithmic scale (10-base). 

t DI are shown in Fig. El The ratio t DI /t TI is shown in Fig. ^2 

All work was done in Matlab 7.0.1, on a 1.66 GHz PC (Intel Core Duo), with 512 MB 
Memory. 

V. DISCUSSION: PROBLEMS AND POSSIBLE SOLUTIONS 

Although the planar TI-FFT algorithm has so many advantages given above, some prob- 
lems do exist in the practical applications. 



1. Complicate geometry 

As an example, consider surface S shown in Fig. where the surface itself is not a 
quasi-planar surface and the direct implementation of the planar TI-FFT algorithm requires 
a large number of FFT operations, which can be seen from the spatial reference planes 
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FIG. 11: The efficiency of the planar TI-FFT algorithm: the ratio of t DI /t TI for J\f x = N y = 
(128,256,512,1024). 

with a spatial slicing spacing 5 Z . The problem can be solved by dividing surface 5* into two 
surface patches A Si and AS2, which can be considered as quasi-planar surfaces and the 
planar TI-FFT can be used on them independently, with coordinate systems selected based 
on the spatial reference planes. At the extreme limit where surface patches A Si and AS2 
are planes, the number of FFT operations reduces to J\f FFT = 2. 

2. Observation points not on the computational grid 

It is well-known that the FFT requires an even grid spacing (but 5 X and 5 y need not to be 
equal), which raises the question of how to calculate the electric field at points that are not 
exactly on the computational grid, e.g., the red filled circles in Fig. One solution for this 
problem is to zero-pad the computational grid in the spectral domain, which corresponds to 
the interpolation of the computational grid in the spatial domain, as shown in Fig. in 
the above example, it has been assumed that the observation points are evenly distributed 
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FIG. 12: An example of complicate surface S that can be divided into two quasi-planar surface 
patches A Si and AS2. The computations of each surface patch is done in its corresponding 
coordinate system whose z-coordinate is perpendicular to the slicing spatial reference planes. 

and the interpolation results are exact provided that the sampling rate is above the Nyquist 
rate |7j. For complicate observation point configurations (e.g., unevenly distributed points), 
the approximate techniques like the Gauss's forward/backward interpolations can be used. 

3. The translation in spatial domain 

In the real situation, the source field surface and the observation surface are separate 
far away from each other (see Fig. H5jl . It is not practical nor necessary to use a large 
computational grid that covers both the source field surface and the observation surface. 
This kind of problem can be solved by using two computational grids, one for the source 
field surface and the other for the observation surface, with the same grid spacings (5 X , S y ). 
Then the translation of the observation coordinate system in the spatial domain, which is 
denoted as (x , yo), corresponds to the phase shift in the spectral domain. Suppose the 
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FIG. 13: The problem of computation of electromagnetic field on the observation points that are 
not on the computational grid (4x4), which are denoted as red filled circles in the spatial domain 
(assume that they are evenly distributed). (6k x , 5k y ) are grid spacings in the spectral domain. (5 X , 
8y) are grid spacings in the space domain. 
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FIG. 14: The zero-padding in the spectral domain (4x4 — > 8x8) corresponding to the interpolation 
in the spatial domain (4 x 4 — > 8 x 8). (5k x , Sk y ) are still the same after zero-padding. But grid 
spacings in the spatial domain become (5 x /2, 5 y /2) after interpolation. 



electric field in the source coordinate system is expressed as E(*< - W - m ), according 
to the property of the Fourier transform [7], the electric field E(x, y) in the observation 
coordinate system is given as 
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FIG. 15: The translation of the source coordinate system o'(0,0) to the observation coordinate 
system o(xq, yo) in the spatial domain. Both the source and observation coordinate systems 
should have the same grid spacings (S x , 5 y ). 



(70) 

4- Computational redundancy 

In the numerical implementation of the planar TI-FFT algorithm, the spatial domain or 
the spectral domain are divided into many small subdomains where the FFT can be used to 
interpolate the electromagnetic field (see Fig. H an d Fig. |SJ). However, the FFT operation is 
done on the whole spatial or spectral domain even though the interpolation is only necessary 
on the relatively small subdomain, which causes the computational redundancy in the planar 
TI-FFT algorithm. Fortunately, the computational redundancy is small for a quasi-planar 
surface and a narrow-band beam. 



E(x )2 /) = IFT 2 



-jkxXo e -jk y y prj> 



E(x' 



xo,y - yo) 
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VI. CONCLUSION 



In this article, the optimized planar TI-FFT algorithm for the computation of electro- 
magnetic wave propagation has been introduced for the narrow-band beam and the quasi- 
planar geometry. Two types of TI-FFT algorithm are available, i.e., the spatial TI-FFT 
and the spectral TI-FFT. The former is for computation of electromagnetic wave on the 
quasi-planar surface and the latter is for computation of the 2D Fourier spectrum of the 
electromagnetic wave. The optimized order of Taylor series used in the planar TI-FFT al- 
gorithm is found to be closely related to the algorithm's computational accuracy 7 TI , which 
is given as Af° pt ~ — ln7 TI and the optimized spatial slicing spacing between two adjacent 
spatial reference planes only depends on the characteristic wavelength A c of the electro- 
magnetic wave, which is 8° pt ~ J7^c- The optimized computational complexity is given as 
Af° p W° pt O (N log 2 N) for an N = Af x x Af y computational grid. The planar TI-FFT algo- 
rithm allows a low sampling rate (large sampling spacing) required by the sampling theorem. 
Also, the algorithm doesn't have the problem of singularity. The planar TI-FFT algorithm 
has applications in near-field and far-field computations, beam-shaping mirror system de- 
signs, diffraction and scattering phenomena, millimeter wave propagation, and microwave 
imaging in the half-space scenario. 
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